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Abstract 

To ease analysis and simulation we make low-dimensional models 
of complicated dynamical systems. Centre manifold theory provides 
a systematic basis for the reduction of dimensionality from some de- 
tailed dynamical prescription down to a relatively simple model. An 
initial condition for the detailed dynamics also has to be projected onto 
the low-dimensional model, but has received scant attention. Herein, 
based upon the reduction algorithm in [27], I develop a straightfor- 



ward algorithm for the computer algebra derivation of this projection. 
The method is systematic and is based upon the geometric picture un- 
derlying centre manifold theory. The method is applied to examples 
of a pitchfork and a Hopf bifurcation. There is a close relationship 
between this projection of initial conditions and the correct projec- 
tion of forcing onto a model. I reaffirm this connection and show how 
the effects of forcing, both interior and from the boundary, should be 
properly included in a dynamical model. 
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1 Introduction 

Ordinary differential equations or partial differential equations are used to 
describe dynamics in the physical world. But many modes of behaviour are 
of little physical interest in particular applications. The essential dynamical 
behaviour of the system is then determined by the evolution of a subset of 
the possible modes; for example, the flight of a ball is dictated solely by its 
velocity and spin, and hardly at all by any internal visco-elastic vibrations. 
Typically we say that the state of the system u is some function, 

u = v(s), (1) 

of the few interesting modes s. Then a model of the dynamics, the rigid- 
body dynamics of a ball for example, is governed by the low-dimensional 
dynamical system 

s = g{a), (2) 

where the overdot denotes d/dt. Rapid damping typically characterises the 
modes that need to be eliminated from consideration. When many modes are 
heavily damped, trajectories are rapidly attracted to some low-dimensional 
subset of the state space, parameterized as in ([!]), a so-called invariant mani- 



fold f33], § 1 . 1 C] . This geometric picture of exponential collapse to a smooth in- 



variant manifold is at the heart of the application of centre manifolds 
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to the rational construction of low-dimensional models (0) by the elimination 
of physically uninteresting fast modes of behaviour. 

The same concepts also lie behind other recent innovations in forming 
low-dimensional models of dynamics: inertial manifolds by Temam |3(J and 
others; the use of centre-unstable manifolds [0, |U, more general invariant 
manifolds 0, |24|, ^ and the so-called nonlinear Galerkin method [19, 

m 



18, 12 



In the earlier paper [^7| I reported an iterative method, based directly 
upon the residuals of the governing differential equations, for the construction 
of such low-dimensional, dynamical models. The evaluation of the residuals is 
a routine algebraic task which is easily programmed into computer algebra; 
programmed without having to become involved in all the messy details 
of asymptotic expansions. Two examples were discussed in [p7[] : a model 
displaying the pitchfork bifurcation of a modified Burger's equation; and the 
long-wave lubrication dynamics of a thin film of fluid. 

However, it is not sufficient to just model the dynamics. To form a 
complete problem, an initial condition, say s(0) = s , has to be specified for 
the model (0). For very low-dimensional models with simple attractors there 
is perhaps little motivation, which may explain the lack of attention given 
to this issue. However, for models of higher dimension, initial conditions 
have long-lasting effects and need to be modelled correctly. This is seen in 



examples such as: the approach to limit cycles [34, 15]; the quasi-geostrophic 
IT], e.g.]; long-wave models of fluid films |28], e.g.] and of 
and the concept of "initial slip" in some 



approximation 
dispersion in channels 
disciplines @, |]| [13 
In 



20 



e.g-J 



|22| , §2] I introduced a simple dynamical system to illustrate the correct 
treatment of initial conditions when forming a low- dimensional model. I now 
summarise that example. Consider the 2-D dynamical system 



x 



y 



-xy. 



X 



(3) 
(4) 



Linearly, y decays exponentially quickly to leave x as the dominant mode 
appearing in the long-term evolution. Nonlinear theory, see Carr [0] for 
example, asserts that the long-term evolution actually takes place on the 
parabolic centre manifold, Ai c , described parametrically as 



x 



y 



(5) 



Exponentially quickly solutions of (|3|-f|) approach solutions on Ai c whose 
evolution is described by the low- dimensional model 
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The question is: given some initial state for the detailed dynamical system (|3|- 
fD, say (xo,y ) not on Ai c , what initial value of s, say s , is appropriate to 
use with the model (0) in order to guarantee fidelity between predictions of 



the model and the original? In |22|, §2.1] I showed that this is achieved by 
setting 



s = x - x (y - x ) + O (e J , (7) 

where e measures the distance from the given initial condition to the centre 
manifold. Alternatively, in a form we use in later generality, this projection 
from the given (x ,y ) to the correct (s , Sq) on M. c is orthogonal to the 
vector z(s) ~ (1, — s). Only then does the predictions of the model accurately 
describe the long-term evolution of the original. 

How in general do we find initial conditions to use with a model to make 
correct long-term forecasts? The Relevance Theorem of centre manifold the- 
ory 0, |27j assures us that there is indeed a particular solution of the low- 



dimensional model on Ai c which is approached exponentially by every tra- 
jectory of the full dynamical system. As developed in p2| , the geometric 
picture of evolution near the centre manifold suggests a method of analysis. 
The algebra is based upon how trajectories near to the centre manifold evolve, 
identifying which ones approach each other exponentially quickly and thus 
have the same long-term evolution. For example, in the long-term disper- 



sion down a channel EQ, 32 1 we can discern the difference between dumping 



contaminant into the slow moving flow near the bank and into the fast core 



flow. Normal form transformations [ 10 ] also support this projection of initial 



conditions onto the centre manifold. In this paper I simplify the main results 
of [^2|, pp65-6], §§2.1|, and in §§|Q| show how to solve the new equations using 



an iterative scheme based on that developed in [|27| to derive the dynamical 
model. 

However, the centre manifold based analysis reported to date has always 
dealt with the issue of initial conditions for models based upon slow modes 
with near zero growth rate. In general, a centre manifold model will involve 
the oscillating dynamics of a Hopf bifurcation or of the loss of stability to 
travelling waves. In §^ we also consider what extensions need to be made to 
the analysis and the iterative algorithm to provide correct initial conditions 
to models of such oscillations. 

Lastly, one attribute of basic centre manifold theory is that it applies to 
autonomous dynamical systems. But the modelling of dynamics with time- 
dependent forcing is of considerable interest. In the presence of forcing a 
system is pushed away from the centre manifold as quantified by Cox & I 



[ 10 ] . Thus there is a close connection between the geometric projection of ini- 



tial conditions and the appropriate projection of forcing onto the model G2 
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§7]. There is also many interesting issues in the modelling of noisy dynam- 
ical systems, as expressed by stochastic differential equations 0. In §§f2~^ 
I reaffirm the connection and show that a new normalisation of the initial 
condition projection onto the centre manifold makes the projection of forcing 
significantly simpler. These are applied in §§ |3.3| to the modified Burger's dy- 
namics to show how both interior and boundary forcing are treated. These 
methods should apply to interesting questions such as: what influence may 
turbulence have on dispersion? and how does substrate roughness affect the 
flow of thin films? 



2 Geometric basis of initial condition projec- 
tion 

The aim of this section is to introduce some of the concepts and details of 
modelling initial conditions of dynamical systems. The presentation improves 



on earlier work reported elsewhere, predominantly in |22|, |8], [|, and is adapted 
to a simple computational approach. 

Consider a general autonomous dynamical system written in a form rel- 
ative to some fixed point (taken to be the origin without loss of generality): 

u = £u + f{u,e), (8) 

where u(t) is the state vector, which may be finite or infinite dimensional, 
£ is the linear operator, and / denotes all the nonlinear terms in the dy- 
namical prescription with possible small parameters e. We suppose that C 
has m eigenvalues with zero real-part and the remaining eigenvalues have 
strictly negative real-part. The linear dynamics of it = Cu then collapse 
exponentially quickly onto the centre subspace £ c spanned by the eigenvec- 
tors and possibly generalised eigenvectors of C, namely e°. Nonlinear theory 
HH correspondingly asserts that the exists an m-dimensional, exponentially 
attractive centre manifold A4 C which has S c as its tangent at the origin. The 
centre manifold is parameterized by any convenient set of m parameters, say 
s. Thus the low-dimensional model of ([8]) is that 

u = v(s,e), such that s = Qs + g(s, e) , (9) 

where v describes the shape of Ai c , Q is a linear operator with eigenvalues all 
of real-part zero, and g is strictly nonlinear. Theory also asserts this model 
is valid exponentially quickly. 



In the earlier paper [27] I derived a robust and straightforward iterative 



algorithm to construct approximations to the model functions v and g. I 
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take this work as read. We now move on to develop the correct projection 
of initial conditions for (|8]), and then how to project a forcing superimposed 
on the autonomous dynamical system. 

For simplicity I first assume that the critical eigenvalues of £ are all zero. 
The case of non-zero, pure imaginary eigenvalues has extra complicating 
details and is addressed in §f| where we investigate the specific example of a 
Hopf bifurcation. 

2.1 Evolution near the centre manifold 

As in §5.1], consider the evolution of a point on the centre manifold, Ai c , 
and the evolution of any neighbouring point. Let n denote the small vector 
joining these points, then, under the flow of the dynamics of (@), n satisfies 
the linear equation 

dft 

— = Jn, where J = C+M, (10) 

is the Jacobian of (|]) evaluated on M. c and where M = df / du is the Jacobian 
of the nonlinear and parameter dependent terms. Note that M and hence J 
are functions of position s on the centre manifold. It is useful to imagine any 
given n to be a function of position on Ai c rather than time; in this case we 
deduce by the chain rule that 

using the summation convention hereafter. The projection of initial condi- 
tions onto A4 C is done along what have been termed "isochronic manifolds" , 
denoted herein by X. Such a projection of initial conditions is also supported 
by normal form transformations as discussed Cox & Roberts ||10|| . The most 
well known isochronic manifold is simply the stable manifold of the origin, 
M. s : since the origin is an equilibrium in the model, that trajectories starting 
on M. s exponentially quickly approach the origin is precisely the requirement 
for the isochronic manifold of the origin's fixed dynamics. 

Under the evolution in the neighbourhood of A4 C the displacement vector 
n will do either of two things: the tip of n off the manifold will approach 
A4 C exponentially quickly but in general it will have slipped from the base 
point on M. c and thus n will exponentially quickly become tangent to A4 C ; 
alternatively, if n is aligned just right then the tip will approach the base 
while remaining transverse to M. c . It is this latter case that is of interest 
as an initial condition of the full dynamics at the tip of n will result in a 
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long-term evolution that is indistinguishable from that of the evolution of 
the base point. Thus the base point forms the appropriate initial condition 
for the low- dimensional dynamics on Ai c . The set of all such tip points that 
exponentially approach Ai c , while n remains transverse, forms the isochronic 
manifold X for any specified initial condition on M. c . 

To describe the projection of given initial conditions tto of the full dynam- 
ics onto A4 C , we use the normal vectors to the isochronic manifolds. More 
specifically we use the normal vectors of X at M. c . Let n a = u — v(s , e) be 
a family of small vectors which span the tangent space of X. Then in terms 
of an inner product ( , ), we seek m linearly independent vectors, say as 
a function of position s, such that 

(r t ,n a ) = 0. (11) 

It is much easier to find the m vectors than the possible infinity of vectors 
n a . Taking d/dt of ( [LTD an d using ( p~Q|) we deduce 

Vri = where V = 4- + J* , (12) 

at 

and denotes the adjoint in the specified inner product. This equation 
describes how the normal vectors vary over _M C - Call T> the dual. 



We need to solve (12). But in general will inconveniently vary quickly 
in magnitude and possibly direction over A4 C - Whereas all we are actually 
interested in is the space spanned by r i: namely the tangent space to X at 
A4 C . Instead we seek an equivalent basis for X, one which varies relatively 
slowly over A4 C , by the invertible transformation = QijZj for some basis 
vectors Zj also a function of position on M. c . It is from here that we depart 
significantly from the analysis reported in my earlier work |22|, §5]. For a 
reason that becomes apparent in the next subsection, we seek the particular 
basis such that 

dv 

(z i ,e j )=5 ij where ej(s,e) = — (13) 

OSj 

are local tangent vectors to Ai c based upon the parameterization of A4 C ; 
typically e,j — > e° as (s, e) — > 0. Substituting r. ; = QijZj into ( |T2"D then leads 
to 



~dT" J 



Zj + QijVzj = 



+ QijCDzj, e k ) = upon taking ( , e k ) 
—Qij(Vzj, e k )z k + QijVzj = putting dQij/dt in the first. 
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But Qij is an invertible matrix and thus we must solve 

Vz j -(Vz j ,e k )z k = 0, (14) 

in conjunction with the orthonormality condition (|T^) in order to find the 
basis vectors Zj for the isochronic manifolds X. 

Equation fll4]) has a reasonable interpretation. The second term just 
projects the residual of the dual flTjp onto X, and hence (|14|) requires all 



other components of the residual to be zero. Thus (|14D requires that the 
basis indeed twists with X, but places no restraint on how the basis spans 
X; the orthonormality condition (|H|) closes the problem to give a unique 
solution. 

Once the basis vectors Zj are found, we then solve 

(z j (s o ,e),u o -v(s o ,e)) = for all j (15) 

to determine the projection of a given initial state Uq onto an initial state Sq 
for the model (|9]). This projection is linear in distance away from the centre 
manifold M. c . There will be errors quadratic in the distance. However, 
in many applications the stable manifold M. s is precisely the linear stable 
subspace £ s ; hence at least near the origin we may expect that a linear 
projection onto M. c will be quite good. 

2.2 Model forcing by the same projection 

In this subsection we consider the dynamical system @ with a small, of 
O (e), forcing superimposed. Namely we analyse briefly 

ii = Cu + f(u, e) + sp(u, t, e) , (16) 

for small forcing ep. The forcing could be deterministic, as investigated by 
Cox & I [^, |8|, or stochastic as examined by Chao & I 0. Our aim is 



to transform the forcing of the detailed system (|T6|) into a corresponding 
forcing of the model ([5|). That is, we seek the forced centre manifold and the 
evolution thereon in the form 

u = v(s,e) + ew(s,t,e) + O (e 2 ^J , 
s.t. s = gs + g( y s,e)+eq( y s,t,e) + 0{e 2 ) , (17) 

where w describes the displacement of Ai c and q is our main interest as it 
describes the correct forcing to be used in the model. 
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That the projection of a forcing onto a model is nontrivial was shown in 
22 , §7.1]. There a steady forcing of — e in the y equation (Q) of the simple 



system (H-H) results in the model 



s = — s 3 + es , 

which exhibits the destabilization of the origin in favour of either of two 
fixed points at s — rt-y/e. This result is remarkable in that the response is 
comparatively large, of O (\fe), when the original forcing is normal to S c and 
so usually would be neglected by heuristic arguments. 

In P2L §7] I argued that the projection of initial conditions could be 
used to deduce how to model forcing. This close connection between the 
two processes was supported by further work by Cox & I Here I briefly 
reaffirm the connection and show why the orthonormality condition ([13]) is 
desirable. 

To find an equation for q, simply substitute flTTD into the original sys- 
tem (16) and group all terms linear in e to deduce 

dw „ . , 

— + £q = Jw + p, (18) 

where S = [ej] is the matrix of tangent vectors, and here 

d- d- d- . 

dt = m + dS {Gs + 9h 

due to the direct dependency upon time introduced by the forcing p(u,t). 
Taking (zj, ) of this equation and using the adjoint properties we must have 

j t (Zi,w) - + ( z h e j)lj = {J ] z h w) + (zi,p) . (19) 

Using the orthonormality QI3]), (Zi,ej)qj = qi. Without loss of generality, 
choose the parameterization of positions near Ai c so that 

(z u w) = 0. (20) 

Indeed, Cox &; I || showed that this is the only choice for w that removes 
clumsy history dependent integrals from the forcing q of the model. Thus fll9]) 
becomes 

dz- 

Qi = ( z i,P) + (-JT + ^ Zi ^ w ) ■ 
at 

The last term involves T>Zi which, by the projected dual (0), must lie in 
the space spanned by the z^s, is thus orthogonal to w, and so the last term 
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vanishes. Hence the appropriate linear approximation to the forcing of the 
model is simply the projection 

Qi = (Zi,p), (21) 

for any given forcing p in terms of the vectors Zi determined for the projection 
of initial conditions. 



3 Initial conditions — a pitchfork bifurcation 



In [57] I described a simple and robust iterative scheme for the computer al- 
gebra derivation of the dynamical model @ from the detailed system In 
this section I describe how to extend the iterative scheme to derive the projec- 
tion of initial conditions. This scheme is also eminently suitable for computer 
algebra and I illustrate its application by using the iteration to determine ini- 
tial conditions for the relatively simple pitchfork bifurcation dynamics in a 
specific infinite dimensional dynamical system. Here we consider the sim- 
pler case of centre manifold models formed when the critical eigenvalues are 
precisely zero, rather than the more complicated case of no n- zero imaginary 
part considered in the next section. 

A summary of the iteration for the centre manifold model is as follows. 
Suppose we know an approximation to M. c and the evolution thereon, namely 

wsiti(s, e) s.t. s Qs + g(s, e) . 

For example, usually we start the iteration with the linear approximation: 
v = S°s and g = 0. Then we seek an improved description, that 

u m v + v' s.t. s^Qs + g + g', 

where primes indicate small correction terms to be determined. Substitut- 
ing into the governing differential equation (|8]), neglecting products of small 
quantities, and approximating coefficients of primed quantities by their ze- 
roth order approximation, we deduce that the corrections satisfy 

— (Gs + g)+ £°g' + °—gs = Cv + Cv' + f(v, e) . 

OS OS 

It is not obvious, but provided the definition of amplitudes are arranged so 
that Q is in Jordan form, we may significantly simplify the algorithm by also 
neglecting the term ^nfQs. (It is often physically appealing to use the Jordan 
form because, for example, the two amplitudes involved often represent the 
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"position" and "momentum" of a specific mode.) Thus, rearranging and 
recognising that 

dv dv 

dS {gs + 9)= ti 

by the chain rule for the current approximation, we solve 

Cv'-£°g' = — -Cv-f(v,e). (22) 
dt 

Recognise that the right-hand side is simply the residual of the governing 
equation (H) evaluated at the current approximation. Thus at any itera- 
tion we just deal with physically meaningful expressions; all the complicated 
rearrangements of asymptotic expansions as needed by earlier methods are 
absent. All the messy algebra in the repeated evaluation of the residuals may 
be left to the computer to perform — such mindless repetition is ideal for a 
computer — whereas all a human need concern themselves with is setting up 
the typical solution of 

Cv' — £°g' = residual , 

and not at all with the detailed algebraic machinations of asymptotic expan- 
sions. 

Consider the following variation to Burger's equation featuring growth, 
(1 + e)u, nonlinearity, uu x , and dissipation, u xx . 

du , . du d 2 u , , , 

— = (l + e ) u + u _ + _ , u(0,f) =«(*,*) = 0, (23) 

for some function u(x,t). View this as an infinite dimensional dynamical 
system, the state space being the set of all functions u(x) on [0, tt}. The above 
iteration scheme may be employed to find the centre manifold dynamics near 
the bifurcation that takes place as e crosses zero. This application is described 
in detail in [27 : lines 1-29 of the reduce^ computer algebra program listed 
in §§ |3.2| tell us that the centre manifold is 

u = a sin(x) + 1 - t ^ a 2 sin(2x) + 1 - 7 ^ 12 a 3 sin(3x) + O (e 2 , a 4 ) , (24) 

when parameterized by a, the amplitude of the sin(x) component of u. Upon 
this centre manifold the low-dimensional model is that 



l-e/3 3 
ea a 



+ 0(e\a i ). (25) 



12 

This model, for example, predicts the pitchfork bifurcation as e crosses zero. 



1 At the time of writing, information about reduce was available from Anthony 
C. Hcarn, RAND, Santa Monica, CA 90407-2138, USA. E-mail: reduce@rand.org 
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3.1 Iterative algorithm for the projection 

Having outlined the iteration scheme to find the centre manifold model, we 
now turn to implementing a similar but novel iteration scheme to determine 
the vectors Zj that govern the projection of initial conditions. The scheme 
is illustrated by applying it to the pitchfork bifurcation in (|23|) . 

We need to solve ( |i4|) subject to the orthonormality condition (0). The 
method of solution is to iteratively improve an approximation based upon 
the residuals of the equations. We start the iteration with the linear ap- 
proximation that Zj are eigenvectors or generalised eigenvectors, z®, of C\ 
the adjoint linear operator. Suppose that at some later stage we know an 
approximation Zj ~ Zj, we then seek an improved approximation 

Zj « Zj + z'j . (26) 

Firstly, substituting into the orthonormality condition (^3|) gives 

(z'^ej) = Sij - (z h ej) . 

Approximating the coefficient of the primed correction quantity then shows 
that we impose on z\ the requirement that 

(4, e?> = - {z h ej) . (27) 

Secondly, substituting ( pE| ) into the projected dual equation (]TJ]) and drop- 
ping products of correction terms gives 

- Vz'j + (Vzj, e k )z' k + (Vz'j, e k )z k = Vz 3 - (Vz j7 e k )z k . (28) 

Approximating all coefficients of primed quantities, all appearing on the left- 
hand side, by their zeroth-order approximation, we seek a correction such 
that 

-tfz'j + {C)z% e°>4 + (tfz'j, el)zl = V~ Z j - (Vzj, e k )z k . 

The two inner products on the left-hand side vanish as they both may be 
transformed to a form ( , Ce k ) which is zero as e k is a critical eigenvector of 
C. Thus we solve the linear equation 

-C)z'j =V~Zj- (Vzj,e k )z k , (29) 

for the corrections z'-. Thus the corrections are simply driven by the residual 
of the projected dual (|H1 ) evaluated at the current approximation as appears 
on the right-hand side of (p9|). One uncomfortable feature of (^) is that 
sometimes during the course of the iteration there is a component of z k in 
the right-hand side — it should be ignored and ultimately it will vanish as the 
iteration proceeds. 
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3.2 The example pitchfork bifurcation 

A principal reason for adopting this approach is because the iteration is 
simply implemented in computer algebra. I discuss the implementation of 
the iterative algorithm when applied to determining initial conditions for the 
model fl25|) of the modified Burger's equation (p3j). 

Based upon the above derivation, the general outline of the algorithm is: 

1. find the centre manifold and the evolution thereon; 

2. initialisation and linear approximation; 

3. repeat until residuals are small enough; 

(a) compute normality and adjoint residuals, 

(b) compute projected adjoint residual, 

(c) solve for the correction and update approximation. 

In practise, the iteration for the initial condition projection could be inter- 
twined with the iteration for the centre manifold model. However, here we 

keep them separate for clarity. 

Implemented in REDUCE for Burger's equation (|23| ) the algorithm may 
look like the second part of the following. 

1 COMMENT First find pitchfork bifurcation in u_t=(l+eps)u+uu_x+u_xx, 

2 where a(t) measures amplitude of sin(x) component in u(x,t) 

3 ; 

4 on div; off allfac; on revpri; factor sin; 7, improve print appearance 

5 % trigonometry rules OK 

6 let { sin(~x) *cos (~y) => (sin(x+y)+sin(x-y) )/2 

7 , sin(~x) *sin(~y) => (-cos(x+y)+cos(x-y))/2 

8 , sin(~x)~2 => (-cos(2*x)+l)/2>; 

9 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/ 
/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/^ 

10 °/ Define the inverse operator of u+u_xx 

11 operator linv; linear linv; 

12 let linv(sin(~k*x) ,x) => sin(k*x)/(l-k~2) ; 

13 % Using inner product: <u,v>=(2/pi)\int_0~pi u.v dx 

14 operator mmean; linear mmean; 

15 let { mmean (l,x)=>2, mmean (cos (x) ,x)=>0, mmean(cos (~k*x) ,x)=>0 }; 

16 7. 

17 depend a,t; 7o asserts that a depends upon time t 

18 let df (a,t) => g; 7o so da/dt is replaced by current g(a,eps) 

19 7. 

20 u:=a*sin(x); g:=0; 7o initial approximation 

21 7. 

22 7o iterate until PDE is satisfied (to requisite order) 

23 let -[eps~2=0, a~4=0}; % discard high-order terms in a & eps 
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24 repeat begin 



25 write eqn:=df (u,t)-(l+eps)*u-u*df (u,x)-df (u,x,x) ; 

26 gd:=-mmean(sin(x)*eqn,x) ; 

27 write u:=u+linv(eqn+sin(x)*gd,x) ; 

28 write g:=g+gd; 



29 end until eqn=0; 

31 % Second find projection of initial conditions onto centre manifold 

32 % Get tangent vector to centre manifold for normalisation 

33 write es:=df(u,a); 

34 7 Define the adjusted inverse of adjoint operator of u+u_xx 

35 operator lainv; linear lainv; 

36 let {lainv(sin(~k*x) ,x) => sin(k*x) / (l-k~2) , lainv(sin(x) ,x)=>0>; 

37 °/ linear approximation to projection kernel, then iterate 

38 z : =sin(x) ; 

39 repeat begin 



40 write norm:=mmean(z*es,x)-l; 

41 dz:= df (z,t)+(l+eps)*z-u*df (z,x)+df (z,x,2) ; 

42 write eqn:=dz-mmean(dz*es,x)*z; 

43 write z:=z-lainv(eqn,x)-norm*sin(x) ; 



44 end until (eqn=0) and(norm=0) ; 

45 end; 

Observe the how lines 32-44 of this program implements the algorithm 
for the initial condition projection. 

1. £4-29 find the centre manifold and the evolution thereon using the 
iterative algorithm of |2TJ as outlined at the start of §|3|; 

2. £33-38 Initialisation. 

• £33, the local tangent vector to Ai c , namely 

1 3 
e(x) = sin x + -a sin(2a;) + —a 2 + O (a 3 + e 3/2 ) , 

• £35-36 defines lainv, the inverse of the adjoint operator Here, 
under the obvious inner product 

2 r n 

(u, v) = — / uvdx , 

7T JO 

£ is self adjoint so lainv is identical to linv, except that we need 
to neglect any component in z°(x) = e°(x) = sinx as commented 
upon earlier. 

• £38 gives the initial linear approximation to the projection 

z ~ z° = sin x . 
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3. £39-44 perform the iteration until the residuals are negligible according 
to the let statement of £23; 

(a) £40-41 compute normality (|13"D and dual fli~2D residuals for the 
current approximation direct from their equations, 



(b) £42 computes the projected dual residual flTj), 

(c) £43 solves for the correction and updates the approximation to the 
projection vector z(x). 

Running this program shows that 

l + e/3 , N l + 5e/4 , , . 

smx asm(2x) H a sm(3x) 

6 96 




(30) 

The task of finding the correct initial condition for the model ( f25|) is then the 
following. Given an initial condition for the Burger's equation (p3|), namely 
that u = uq(x) at t = 0, we project onto the centre manifold by solving for 
the amplitude ao in the nonlinear equation 

(z(ao, x), Uo(x) — v(ao, x)) — . (31) 

An iterative approach will usually suffice to solve this. Starting with the 
approximation ao = (z ,Uq), successive corrections may be computed as 

a' = (z(d , x),u (x) - v(a , x)) . 

For example, if uq = a sin x for some particular a, then 

a = a + ^a 3 + O (a\ e 2 ) , (32) 

and not simply ao = a as would be implied by a direct application of the 
definition of amplitude a. 

Of course in this particular application the issue of the precisely correct 
initial condition is of little interest because the ultimate fate of the dynamics 
is absorption by a stable fixed point and an incorrect initial condition just 
causes a small error in the timing of the approach. However, in more compli- 
cated dynamical models with non-trivial long-term dynamics, for example in 
chaotic models or in the shear dispersion of contaminant in a pipe or chan- 
nel, errors in the initial condition can cause significant long-term errors in 
the predictions of a model. But before moving on to the analysis of such a 
problem, we investigate in the next subsection the projection of forcing in 
this modified Burger's equation (^3]). 
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3.3 Forcing in the interior and on boundaries 

The forced equation we consider briefly in this subsection is ( |23"D with some 
unspecified small forcing sp(u,x,t), namely 

du du d 2 u 

- = {l + e) u + u - + — + ep{u,x,t). (33) 



Then by the arguments of §2.2j the forcing of the model fl25|) turns it into 

l-e/3 



a = ea 



12 



-a 3 + (z, epiy, x, t)) + O (e 2 , a 4 , e 2 ) . (34) 



For example, a spatially uniform, additive forcing ep(t) induces a multiplica- 
tive forcing as in 




3 , / 4 , 17 + 5e/4 2 \ , 2 4 , 



72tt 



-a 2 \ep + (c z ,a 4 ,e 



Also, as in the example near the start of §2.2j , here a forcing proportional to 
sin (2a;) may destabilize the origin. 

The above results on forcing in the interior of the domain are straightfor- 



ward given the analysis of § |2.2| . A little more subtle is the effects of forcing 
in the boundary conditions. For example, ( |33"D may have forced boundary 
conditions such as 

u(0, t) = ep (t) , u(n, t) = ep n (t) . (35) 

For the moment, assume there is no forcing in the interior — these are the 
only forcing terms. 

One heuristic approach is to turn these boundary conditions into inte- 
rior Dirac delta function forcing of the same problem but with homogeneous 
boundary conditions. In essence, this approach forces an extremely thin 
boundary layer at x = + between the boundary value of u(0) = and an 
interior value of w(0 ++ ) = epo, and similarly near x = tt. We try the forcing 
ep = A5'(x — + ) + B8'(x — 7r~) in ( |3~3D with the homogeneous boundary 
conditions of fl2"3]). Here 5'(x) denotes the derivative of the Dirac delta func- 



tion. Then integrating x times (|33| ) over the extremely small interval [0, ++ ] 
leads to w(0 ++ ) = —A as the effective boundary value of u at x = 0. Similar 
integration near x = tt leads to u(n ) = B as an effective boundary value. 
Thus to match the forced boundary conditions (|35|), we choose 

ep(x, t) = -ep (t)5'(x - + ) + ep v (t)S'(x - tt") . (36) 
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Then the projection of such an "interior" forcing onto the model, equa- 
tion ([34]), leads to 

1 - e/3 3 2 / , 25 2 



ea a H — 1 H a £(»„- + on) 

12 7T V 288 7 yF y ' 

' ' /3 ae(^-p ) + o(e 2 ,a 3 ,£ 2 ) . (37) 



3tt 

Notice that an asymmetry in the boundary forcing, p v ^ po, may destabilize 
the fixed point at the origin. 

Another more systematic approach identifies where inhomogeneous bound- 
ary conditions such as ( |3~5"D affect the earlier analysis. We now do this. Realise 
that with forcing in the boundaries the differential Jacobian term Jw in ([18]) 
then comes with the attached boundary conditions that w(0,t) = po and 
w(n,t) = p n . Consequently, here the integration by parts in going from (|i~8|) 
to ( |T9P introduces extra terms as 

2 

(z,Jw) = (jh,w) [z x w)l . 

7T 

Hence the forcing of the model is not just the projection of the interior 
forcing, (z,p), but instead contains extra terms: 

q=(z,p) + - [z x {0)p {t) - z x (7r)p n (t)} . (38) 
After substituting in the expression (|30|) for z, this agrees precisely with 



4 Correct phase near a Hopf bifurcation 

In the analysis and example of the previous sections, the dynamics on the 
centre manifold have been based on the critical eigenvalues being precisely 0, 
not the more general case where just the real part of the eigenvalues are zero 
but the imaginary part is non-zero. However, because an eigenvalue with a 
non-zero imaginary part is always associated with oscillations, such a case is 
important in practise as it arises in the common transition from steady to 
oscillatory dynamics. Because it is significantly more difficult to determine 
the projection of initial conditions for such oscillatory dynamics, I describe 
a specific implementation in this section. Note that in a Hopf bifurcation, 
as recognised by Winfree |34] and Guckenheimer unless a good initial 



condition is found the phase between the model and the actual dynamics are 
irretrievably different. 
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As an example, consider the dynamics of the following dynamical system 



u 



1 -1 

2 1 
L 2 








•1 



u + 



eui — 2u\Us 



(17, 



(39) 



where e is a control parameter. It is straightforward to discover that when 
e is zero, the critical value, the linear operator has eigenvalues ±i and — 1. 
Thus there exists a centre manifold corresponding to the two eigenvalues ±z. 
The mode with eigenvalue —1 is representative of the many exponentially 
decaying modes we find in real applications. 



4.1 The example centre manifold model 

Our interest herein lies in the provision of correct initial conditions for a 
low-dimensional model, not immediately in the construction of the model. 
Hence, in this subsection I record the principal features of the centre manifold 
model of ([39|) and do not describe its derivation. 

The eigenvectors corresponding to the critical eigenvalues enables us to 
construct the centre eigenspace, S c , the linear approximation to the centre 
manifold. The eigenvectors of A = ±i are (2, —2 ± 2i, —3 ± i) and we use 
the real and imaginary parts of these eigenvectors to span S c . A linear 
approximation to the centre manifold is then 



Sr. 





2 




' " 


where e° = 


-2 


e° = 


2 




-3 


1 



Then in terms of s = (x,y), the linear dynamics on the centre manifold are 
the oscillations described by 



s — Qs , where Q 



-1 

1 



(40) 



In a general problem the centre subspace is described by u — 
the linear evolution equation du/dt = Cu becomes 

k ~ dt ~ j ^ e jyjkSk ■ 



dt 



CjSj. Hence 



Since this holds for all sufficiently small Sk, we deduce the generalised basis 
vectors for the centre subspace satisfy 



Ce* = eSG 



(41) 
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Similarly look at the linear equation for the leading order projection vectors 
z Q j. At leading order the projected dual ([14] ) becomes 

£t z o_ (£ t z )e o )z o = 0) 
after dropping dzj/dt and AH Zj as being small. Then 



{£ z< j> e k) Z k 



{Zj,Ce k )z k by adjoint property 

{z%e\Q lk )zl by© 

8jeGik z °k by orthonormality ( |13|) 



Thus the initial approximation for the projection vectors must satisfy 

Oz® = g jk z° k . 



(42) 



For the nonlinear description, define the "amplitudes" to be precisely 



x 



(z°, u) , y = (z° u) , where 



" 1/2 " 




- 1/2] 





' z y 


1/2 











(43) 



in terms of original variables u. Using the first part of the computer algebra 
program listed in §3~2", a quadratic approximation to the nonlinear shape of 
the centre manifold is 



u 



2x 

-2x + 2y 

-3x + y + fy 2 + ^xy + fx 2 + ey + ex 



+ O (e 3 + x 3 + y 3 ) . (44) 



The lowest order structurally stable model on the centre manifold is the 
following cubic model 



x = —y — 2xy + 6x + ex 

-fxy 2 - fx 2 y - ±fx 3 - 2exy - 2ex 2 + O (e 2 + x 4 + y 4 ) 
y = x + ex + O (e 2 + x 4 + y A ) . 



(45) 



As is generally the case in Hopf bifurcations, with these correct cubic non- 
linearities this model is usefully predictive. Numerical simulations show the 
birth of a limit cycle as e crosses through zero. Here we have systematically 
reduced the dynamics by a small step down from 3-D to a 2-D model. In seri- 
ous applications to very high dimensional dynamics, we would have reduced 
enormously the dimensionality of the dynamics. 
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4.2 Projection onto oscillating dynamics 



In this subsection I first generalise the iterative algorithm presented in §371 
to the case of oscillatory dynamics on the centre manifold. To illustrate the 
algorithm and the typical results I then apply it to the simple dynamical 
system (]39f) . 

Restart the generic nonlinear analysis from (p^). The right-hand side is 
just the residual of the projected dual equation (|14|); it stays the same. The 
operator on the left-hand side has to be simplified, but not as drastically as 
in 



The third term {T^z 1 -, e k )z k only changes that part of the left-hand side 
in the space spanned by {z k }. But we do not solve the dual (|TJ) in 



this space, hence the projection seen in ([14]). Thus this term is safely 
ignored. 

The inner product in the second term of fl28[ ) simplifies under approx- 
imation as follows. 

(VZj,e k ) « (T>z j} e° k ) as e k « e° 

~ (£*z°,e°) neglecting small terms 

= (Gjtz%, el) using (||) 

= G jk as {z% e° k ) = 5 ik . 



The first term, the dual operator, approximates to the homological 
operator 



dz' 

{GkiSe + g^gf + Jiz'j « 



dz' 

GkiSi-7^- + tfz'j 



Thus in general we solve 

dz'- 
r> 3 
— ykeSiT^ — ' 

OSu 



C)z'; + Q jk z' k = Vzj - (DZj, e k )z h 



(46) 



for the corrections to the projection vectors. 

The complicating feature of (HST) is that it is a coupled set of equations 



for the correction vectors z'j — coupled through Gjk z 'k — an< ^ ^ ^ s ^° ^ e 
solved in the space of multinomials in the amplitudes s — because of the 
structure of the homological operator. As noted earlier, in the iteration we 
may ignore components in z^ but we must enforce orthogonality via the 



iterative correction (p7[). Other than these complications the outline of the 
iterative algorithm is the same as in the previous section. 
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For the specific dynamical system (|3^) a reduce computer algebra pro- 
gram follows. Observe that the first part of the program computes the de- 
scription of the centre manifold model. It is the second part that determines 

the projection. 

1 COMMENT Use iteration to form the centre manifold model of an 

2 elementary Hopf bifurcation problem. Then compute the projection 

3 of initial conditions onto the centre manifold by iteration. 

4 ; 

5 X formating for printed output 

6 on div; off allfac; on revpri; factor del; 

7 procedure dfv(v,t); X patch differentiation of vectors 

8 mat((df (v(l,l) ,t)) ,(df (v(2,l),t)) ,(df (v(3,l) ,t))) ; 

9 procedure dfm(z,t); X patch differentiation of matrices 

10 mat((df (z(l,l) ,t) ,df (z(l,2) ,t)) 

11 , (df (z(2,l) ,t) ,df (z(2,2) ,t)) 

12 , (df (z(3,l) ,t) ,df (z(3,2) ,t))) ; 

14 X PART 1: find the centre manifold 

15 X matrix of linear terms and centre eigen-vectors 

16 ll:=mat((-l+del*eps,-l,0) , (2 , 1 ,0) , (1 ,2 , -1) ) ; 

17 ex0:=mat((2) , (-2) , (-3)) ; 

18 ey0:=mat((0),(2),(l)); 

19 X 

20 X only retain terms up to order o in x, y & eps 

21 X use del to count the number of x, y & eps factors in a term 

22 o:=3; let del~4=>0; 

23 X 

24 X linear solution 

25 u:=del*(x*exO+y*eyO) ; 

26 g:=mat((-y) , (x)) ; 

27 depend x,t; let df(x,t) => g(l,l); 

28 depend y,t; let df(y,t) => g(2,l); 

29 X 

30 X iteration 

31 X set arbitrary multinomial & its coefficients for later use 

32 operator c; 

33 hd:=for m:=0:o sum for n:=0:o-m sum c (m,n) *x~m*y~n$ 

34 clist:={>$ 

35 for m:=0:o do for n:=0:o-m do clist:=c(m,n) .clist$ 

36 repeat begin 

37 X compute residual of DDEs 

38 ru:=dfv(u,t) ; ru: =-ru+ll*u+mat ( (-2) , (2) , (0) ) *u(l , 1) *u(3, 1) 

39 +mat((0),(0),(l))*u(2,l)~2; 

40 X solve first two components for evolution g 

41 gd:=mat((ru(l,D), (ru(l,l)+ru(2,l)))/2/del; 

42 write g:=g+gd; 

43 X form and solve homological equation 

44 eqn : =hd-y*df (hd , x) +x*df (hd , y) - (ru (3 , 1 ) - (-3*gd (1,1) +gd (2 , 1 ) ) *del) ; 
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45 elist:={}; for m:=0:o do for n:=0:o-m do 

46 elist :=coef fn(coeffn(eqn,x,m) ,y,n) .elist; 

47 write u:=u+mat((0) , (0) , (sub(solve(elist, elist) ,hd))) ; 

48 end until (ru=mat ( (0) , (0) , (0) ) ) ; 

50 % PART 2: compute the vectors zx & zy to give projection onto CM 

51 % adjoint jacobian and tangent vectors 

52 ja:=tp(ll+mat((-2*u(3,l) ,0 ,-2*u(l , 1) ) 

53 , (2*u(3,l) ,0, 2*u(l,D) , (0,2*u(2,l) ,0)))$ 

54 z0:=mat( (1/2, 1/2), (0,1/2), (0,0)); 

55 es:=mat((df (u(l,l) ,x) ,df (u(l,l) ,y)) 

56 , (df (u(2,l) ,x) ,df (u(2,l) ,y)) 

57 , (df (u(3,l) ,x) ,df (u(3,l) ,y)))/del$ 

58 % truncate approx & form arbitrary multinomials & coefficients 

59 o:=2; let del~3=>0; 

60 zxd:=for m:=0:o sum for n:=0:o-m sum c(m,n, l)*x~m*y~n$ 

61 zyd:=for m:=0:o sum for n:=0:o-m sum c (m,n, 2) *x~m*y~n$ 

62 elist :={}$ 

63 for m:=0:o do for n:=0:o-m do elist : =c(m,n, 1) . (c (m,n, 2) . elist) $ 

64 °/ initial linear approx, then iterate 

65 z:=z0$ 

66 repeat begin 

67 °/ compute residuals of the projected dual and orthogonality 

68 dz:=dfm(z,t) ; dz:=dz+ja*z; 

69 rdz : =dz-z* (tp (es) *dz) ; 

70 rze:=mat((l,0) , (0 , 1) ) -tp(z) *es ; 

71 °/ form and solve homological equation from 3rd component 

72 eqx:= zxd+y*df (zxd,x) -x*df (zxd,y) -zyd -rdz(3,l); 

73 eqy:= zyd+y*df (zyd,x) -x*df (zyd,y)+zxd -rdz(3,2); 

74 elist:={}$ for m:=0:o do for n:=0:o-m do 

75 elist :=coef fn(coeffn(eqx,x,m) ,y,n) . 

76 (coef fn(coef fn(eqy ,x,m) ,y,n) .elist) ; 

77 csoln:=solve(elist, elist) ; 

78 °/ update approx 

79 write z:=z+z0*tp(rze) 

80 +mat((l) , (-1/2) , (l))*mat((sub(csoln,zxd) , sub(csoln,zyd) ) ) ; 

81 end until (rdz=mat ( (0 ,0) , (0, 0) , (0 ,0) ) ) and(rze=mat ( (0 ,0) , (0 ,0) ) ) ; 

82 end; 

Observe how the second part of this program implements the algorithm, 
in lines 50-81. 

1. £14-48 compute the centre manifold ( fHf ) and the model evolution (|45|) 



using the iterative algorithm of [27|. Note that reduce does not im- 



plement matrices well and so £7-12 define two procedures to help. 
2. £51-65 is initialisation. 

• £52-53 computes the adjoint operator on M. c - 

• £54 sets the linear projection vectors z® into a 3 x 2 matrix. 
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• £55-57 computes the tangent vectors ej of M. c . 

• £59 says to truncate the computations to have errors O (<5 3 ) as 
these should be accurate enough and will speed computation. 

• £60-63 define general multinomial expressions and a list of their 
coefficients for later use in solving the homological equation via 
the method of undetermined coefficients. 

• £65 sets Zj to its initial approximation. 

3. £66-81 perform the iterations until the residuals of the desired equa- 
tions are zero to the order of error specified. 



(b) 



£68-70 computes the residuals of the projected dual equation flT4] ) 
and the orthonormality condition (|13"D, assigned to rdz and rze 
respectively. 

£72-77 solves for the undetermined coefficients of the corrections 



using the third component of equation (46j). The first term on the 



right-hand side of £70, for example, comes from —£)z'-\ the second 

dZ' 

and third terms represent —Qkl s t~of-\ while the fourth term comes 
from Gjkz'k'i an d the last term is from the residual. 

(c) £79-80 updates the approximation to better satisfy the orthonor- 
mality and the projected dual. 

The output of this computer program indicates that projections from 
initial conditions off M. c onto M. c are to be orthogonal to 
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(47) 



(48) 



For example, consider the dynamics from the initial condition uq = 
(0.022, 0, 0.073) and the corresponding evolution on M. c . According to either 
the leading order projection or the definition of the amplitudes x and y, equa- 
tion ([43|), Uq corresponds to the point on Ai c with parameters x = y = 0.011, 
namely Uq = (0.022,0,-0.019). However, at this point on Ai c the projec- 
tion of initial conditions should be slightly different; according to the linear 
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modifications in (J7-S|) it should be orthogonal to the columns of 



0.492 
0.004 
-0.008 



0.496 
0.502 
0.004 



Thus, according to (p~5|) , a more appropriate initial condition on A4 C is 
Mq = (0.021,0.001,-0.017). Similarly, the second order corrections in (f47|- 



48]) refine the initial conditions further. The differences here are rather small, 
but the improvement is seen in Figures [l]-^| Figure [l] shows the qualitative 
picture of the trajectory starting at no exponentially quickly approaching 
A4 C , but that the model is only quantitatively predictive if Uo is projected 
correctly. Figure |2] quantifies a comparison between the various orders of ap- 
proximation to the projection and demonstrates that the refined projection 
vectors Zj do perform better. 



Also recall from § |2.2| that even such small corrections to the projection 
as we see here may have a considerably larger influence when projecting an 
applied forcing onto the centre manifold model. 



5 Concluding remarks 

Work in progress will show how to apply the techniques presented here to 
problems of more physical interest. In particular I am examining the issue of 



projecting initial conditions onto the lubrication model |28| of the flow of a 
thin film of fluid over a solid substrate. Although there are no inertia effects 
in the lubrication model of the dynamics, there may well be such effects in 
the provision of initial conditions as the fluid dynamics relaxes to lubrication 
flow. Additionally, the effects of substrate roughness on the flow may be 
determined by projecting the appropriately perturbed substrate boundary 
conditions. 

Finding and coding the adjoint can be a major headache, especially 
for problems such as free-surface fluid dynamics. I have sought approaches 
based directly upon the Jacobian J rather than its adjoint, because then the 
residual driving the iteration could be determined directly and very simply 
from the governing equations. However, so far, the only method I have found 
also involves determining the equivalent of w, introduced in (|i~7D , which is 
considerably more involved. 

Throughout this paper we have addressed the projection of initial condi- 
tions and forcing, linearly correct in distance from the centre manifold or in 
forcing amplitude. If one needs to determine effects nonlinear in distance, 



5 Concluding remarks 



25 




"1 



Figure 1: trajectories of ( |39| ) from the "true" initial condition of Uq (solid), 
and from the various projections onto M. c : linear from z° (dotted); first order 
(dot-dashed); second order (dashed). Observe the approximately exponential 
approach between the "true" and the model trajectories, but that the linearly 
projected trajectory is slightly awry. 
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Figure 2: log of the separation between the trajectories shown in Figure p] 
and the "true" solution: linear from z° (dotted); first order (dot-dashed); 
second order (dashed). Observe the initial exponential approach, but that 
the first order and second order approximations to the initial conditions of 
the model are an order of magnitude better. 
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then a more sophisticated analysis is needed. At this stage the only approach 
I can imagine also involves the considerable labour of finding w. 

Lastly, Muncaster and Cohen |2T], |6|] suggested the construction of the 
low-dimensional manifold of slow, rigid-body dynamics by neglecting rapidly 
oscillating modes. An extremely simple example of the motion of a one-di- 
mensional elastic body is discussed in |25|, §2]. In contrast to the rapid col- 
lapse to the centre manifold, the slow dynamics on the slow manifold form a 
low-dimensional model because they act as a "centre" for the fast oscillations 
of neighbouring trajectories. This principle of neglecting fast oscillations is 
precisely equivalent to the guiding centre principle of Van Kampen |JT| . The 
algebra needed to develop slow manifold models is identical to that pre- 
sented here and in [P^| . The algebra to model initial conditions will also be 
the same. However, slow manifolds are much more delicate. Cox & I |9|, fL0[] 
used normal forms to show that the dynamics on and off the slow manifold 
generally differ by an amount of O {s 2 ), where e measures the amplitude of 
the fast oscillations, it measures the distance off .Mo- That is, generically 
there is some unavoidable slip between a slow model and the fully detailed 
oscillating dynamics. 
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